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ABSTRACT 

We have performed a fully three-dimensional general relativistic magnetohy- 
drodynamic (GRMHD) simulation of jet formation from a thin accretion disk 
around a Schwarzschild black hole with a frcc-falling corona. The initial sim- 
ulation results show that a bipolar jet (velocity ~ 0.3c) is created as shown 
by previous two-dimensional axisymmetric simulations with mirror symmetry at 
the equator. The 3-D simulation ran over one hundred light-crossing time units 
(t's — fs/c where rs = 2GM/c^) which is considerably longer than the previous 
simulations. We show that the jet is initially formed as predicted due in part to 
magnetic pressure from the twisting the initially uniform magnetic field and from 
gas pressure associated with shock formation in the region around r = 3rs. At 
later times, the accretion disk becomes thick and the jet fades resulting in a wind 
that is ejected from the surface of the thickened (torus-like) disk. It should be 
noted that no streaming matter from a donor is included at the outer boundary 
in the simulation (an isolated black hole not binary black hole). The wind flows 
outwards with a wider angle than the initial jet. The widening of the jet is con- 
sistent with the outward moving torsional Alfven waves (TAWs) . This evolution 
of disk-jet couphng suggests that the jet fades with a thickened accretion disk 
due to the lack of streaming material from an accompanying star. 

Subject headings: AGN: accretion disk, jet formation, GRMHD, Schwarzschild 
black hole 



1. Introduction 

Relativistic jets have been observed in active galactic nuclei (AGNs) (e.g. Urry & 
Pavovani 1995; Ferrari 1998; Blandford 2002), in microquasars in our Galaxy (e.g. Mirabel 

& Rodriguez 1999), and it is believed that they originate in the regions near accreting (stellar) 
black holes and neutron stars (Meier, Koide, & Uchida 2001). To investigate the dynamics of 
accretion disks and the associated jet formation, we have performed jet formation simulations 
using a full 3-D GRMHD code. This magnetic-acceleration mechanism has been proposed not 
only for AGN jets, but also for (nonrelativistic) protostellar jets (see Meier, Koide, & Uchida 
(2001)). Kudoh, Matsumoto, & Shibata (1998) found that the terminal velocity of a jet is 
comparable to the rotational velocity of the disk at the foot of the jet in nonrelativistic MHD 
simulations, and that a relativistic jet should be formed near the event horizon. It should be 
noted that the mechanism of jet generation in neutron star binaries may be different from 
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black hole binaries. Recently, an ultra-relativistic outflow was observed from an accreting 
stellar-mass neutron star (Fender et al. 2004). Fender et al. (2004) have concluded that the 
generation of highly relativistic outflows docs not require properties unique to black holes, 
such as an event horizon. It seems plausible that production of ultrarelativistic jets may be 
a generic feature of accretion onto neutron stars such as Cir-X and Sco X-1 at or near to the 
Eddington limit. This result differs from the suggestions that the velocities of jets would be 
limited to the surface escape speed (Livio 1999). 

The fluorescent iron line (K^) seen in the X-ray spectrum of many active galactic nuclei 
(AGNs), and particularly in Seyfert galaxies, has provided us with an extremely useful probe 
of the region close to an accreting black hole. This is an important diagnostic with which 
to study Doppler and gravitational redshifts, thereby providing key information on the local 
kinematics of the material (see the reviews by Fabian et al. 2000 and Reynolds & Nowak 
2003). The simultaneous X-ray /IR flares from a black hole/relativistic jet system are the 
observational evidence revealing a time-dependent interaction between the jet and inner disk 
(e.g. Ueda et al. 2002). In regions where disk-jet interactions take place near a black hole, 
the plasma and magnetic fields interact in a complicated manner; therefore a 3-dimensional 
GRMHD simulation such as we perform here is required to investigate dynamics near the 
horizon and ultimately to couple predicted dynamics to the iron line. 

Koide, Shibata, & Kudoh (1999) have investigated in 2-D, the dynamics of an accretion 
disk initially threaded by a uniform magnetic field in a non-rotating corona (either in a 
state of steady fall or in hydrostatic equilibrium) around a non-rotating (Schwarzschild) 
black hole. The numerical results show that matter in the disk loses angular momentum 
by magnetic braking, then falls into the black hole. The disk falls faster in this simulation 
than in the non-relativistic case because of general-relativistic effects that are important 
below 3rs, where rs = 2GM/c^ is the Schwarzschild radius. A centrifugal barrier at r = 2rs 
strongly decelerates the infalling material. Plasma near the shock at the centrifugal barrier is 
accelerated by the J x B force and forms bipolar jets. Inside this magnetically driven jet, the 
gradient of gas pressure also generates a jet above the shock region {gas-pressure driven jet). 
This two-layered jet structure is formed both in a hydrostatic corona and in a steady-state 
falling corona. Koide et al. (2000) have also developed a new GRMHD code using a Kerr 
geometry and have found that, with a rapidly rotating (a = J/ Jmax = 0.95, J-^ax. = GM"^ /c, 
J and M are angular momentum and mass of black hole) black-hole magnetosphere. the 
maximum velocity of the jet is (0.3 - 0.4) c. Recently the extraction of rotational energy 
from plasma in the ergosphere has been investigated (Koide et al. 2002; Koide 2003). All 
of these GRMHD simulations were made assuming axisymmetry with respect to the 2;-axis 
and mirror symmetry with respect to the plane z — It is known that the axisymmetric 
assumption suppresses azimuthal instabilities. 
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Magnetic fields are expected to play an essential role in the dynamics of most relativistic 
objects such as active galactic nuclei, X-ray binaries, gamma-ray bursts, and quasars. In ac- 
cretion disks, the most important role of the magnetic field is angular momentum transport, 
which twists the field, and results in torsional Alfven waves (TAWs) and associated jet for- 
mation (Koide et al. 1999). Based on the importance of magnctorotational instability (MRI) 
in accretion disks (Balbus & Hawley 1998), the nonlinear growth of these instabilities has 
been investigated by three-dimensional global MHD simulations of black hole accretion disks 
(Hawley 2000; Igumenshchev, Abramowicz, & Narayan 2000; Hawley 2001; Hawley & Kro- 
hk 2001; Machida, Hayashi, & Matsumoto 2000; Hawley, Balbus, & Stone 2001; Armitage, 
Reynolds, k Chiang 2001; Reynolds & Armitage 2001; Hawley k Krolik 2002; Hawley k Bal- 
bus 2002; Krolik k Hawley 2002; Machida k Matsumoto 2003; Igumenshchev, Narayan, k 
Abramowicz 2003) using a pseudo- Newtonian potential (Paczyhsky k Witta 1980). Machida 
k Matsumoto (2003) investigated black hole accretion fiow using three-dimensional global 
resistive MHD simulations. They conclude that the rcconncction in the innermost region 
releases enough energy to explain the origin of X-ray bursts from black hole candidates. 

Recently, more GRMHD simulations have been developed by De Villiers k Hawley 
(2002, 2003a); Gammie, McKinney, k Toth (2003). In particular, Gammie, McKinney, 
k Toth (2003) have employed the numerical scheme HARM (High Accuracy Relativistic 
Magnetohydrodynamics) in which the conservation schemes for GRMHD are used (Koide, 
Shibata, k Kudoh 1998, 1999; Koide et al. 2000). A nonconservative scheme for GRMHD 
following a ZEUS-hke approach has been developed also by De ViUiers k Hawley (2002). 
They also compare the code performance of 1-D and 2-D test cases. De ViUiers k Hawley 
(2003b) have investigated the properties of accretion flows with tori in the Kerr metric using 
3-D GRMHD simulations. The MRI in torus has also been investigated using Schwarzschild 
and Kerr metrics (De Villiers, Hawley, k Krolik 2003; Hirose, et al. 2004, De Villiers, et al. 
2004; Krohk, Hawley k Hirose 2004). 

Different studies have chosen various initial magnetic fields. We have chosen an initial 
magnetic field from the Wald's exact analytic solution for a Kerr black hole immersed in a 
uniform magnetic field (Wald 1974; Punsley 2001; Komissarov 2004). De Villiers, Hawley 
k Krolik (2003) have performed simulations in which the initial conditions for the accretion 
disks and magnetic fields are appropriate for MRI research. They have investigated the 
properties of accretion flows in the Kerr metric including bi-conical outflows in the quarter 
sphere with periodic boundary conditions in the azimuthal direction. These simulations use 
the torus accretion disks as an initial condition. The inner edge of the accretion disk is 
located at r = 15rg, and the inner edge moves into the innermost stable radius within one 
orbit (De Villiers, Hawley, k Krolik 2003). In these simulations, the initial field consists of 
axisymmetric (azimuthal) field loops, laid down along an isodensity surface within the torus. 
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which is different from the initial uniform magnetic field used in simulations by Koide, 
Shibata, & Kudoh (1999) and Nishikawa et al. (2002, 2003). The initial weak loop magnetic 
fields contribute to the generation of the MRI and the associated jet. However, the structure 
of magnetic fields piled up near the black hole in their simulations is different from the 
magnetic fields which are twisted and pinched (not wound) by the free-falling corona and 
accretion disk (Koide, Shibata, & Kudoh 1998, 1999; Koide et al. 2000; Nishikawa et al. 
2002, 2003). The absence of a free-faUing corona (Bondi & Hoyle 1944; Bondi 1952; Michel 
1972) as well as the initial loop magnetic field in (De ViUiers, Hawley & Krohk 2003) results 
in different dynamics of the accretion disk and associated jet formation, compared to the case 
with the free-falling corona and initial uniform magnetic field in Koide, Shibata, & Kudoh 
(1999) and Nishikawa et al. (2002, 2003). 

In this paper we present a fully three-dimensional, GRMHD simulation of jet formation 
with a thin accretion disk. Our motivation for this type of simulation is to determine 
the parameters necessary for relativistic jet formation and the resulting interaction and 
instabilities found between the accretion disk and the black hole. The simulation we present 
here was performed using the same parameters as in Koide, Shibata, & Kudoh (1999) in 
order to determine the physical differences resulting from a full 3-D versus a 2-D simulation 
with axisymmetry and mirror symmetry at the equator. The three-dimensional simulation 
allows us to study the evolution of jet formation because it is run for longer light-crossing 
times than the previous two-dimensional simulations. We find that at the later stages of the 
simulation the accretion disk becomes thick and a wind is formed with a much wider angle 
than the coUimated jet formed at the earlier stage. In Section 2, we describe the governing 
GRMHD equations. The simulation model and numerical method is described in Section 3. 
Our new simulation results are presented in Section 4 followed by a summary and discussion 
in Section 5. 



2. 3-D GRMHD Governing Equations and Numerical Approach 

The governing equations for GRMHD are written in conservation form based on Wein- 
berg (1972) and Misncr, Thorne, & Wheeler (1973). Four-dimensional space-time parame- 
ters are written with Greek indices, while Latin indices are used for three-dimensional space, 
,T° = (ct, x^). The physical parameters include, four-velocity (f/°), proper density (p), proper 
pressure (p). adiabatic index (F), and proper internal energy density (cint = pc^ + p/(F — 1)). 
Relativistic parameters include the Christoffel symbols (F^^) and metric terms {gap)- The 
pressure is related to the other variables by the equation of state, p oc p^. The specific 
relativistic enthalpy is given by, h — Cint + p — pc^ -|- Fp/(F — 1). 
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The conservation equations are derived from the continuity equation, V-y {pU'^) — 0, and 
the conservation of the stress energy tensor, V/.T^'^ = T.^" = + T°"r(^^ + Tf^T^^ = 0. 
From Maxwell's equations, the electromagnetic field-strength tensor is denoted by F^,^ = 
daAf^ — dpAa and the four-current density by J^, where is the four-vector potential. The 
Maxwell equations are written d/^Fap + daFp^ + dpF^a = and VaF"'^ = —J^ . 

For the simulation in the context of this paper, we have assumed infinite electric con- 
ductivity, FajsU" = 0. The stress energy tensor (T"^) for an ideal fluid in an electromagnetic 
field and in a general relativistic environment is, 

T"'' ^(^p+ Pfzri) ^"^^ + P9''^ + - 9'"^^-j^- (1) 

We assume that the off-diagonal spatial elements of the metrics g^i, vanish: gij = (i 7^ j). 
Here Roman indeces run from 1 to 3, the 

5-00 = -hi, Qii = /i- , QiQ = QQi = -hluji/c, (2) 
The line element is written 

3 

^ -a'^i^cdtf + ^{hidx' - c(3'adtf. (3) 

i=l 

Here, the lapse function and shift velocity are defined by a = y^^o + Si=i(^7^)^ /3* ~ 
respectively. 

The 3 + 1 form of equations are employed (Weinberg 1972; Thorne, Price, &; Macdonald 
1986). The 3 + 1 set of equations used in these simulations for the general relativistic con- 
servation laws governing the plasma and Maxwell equations were derived in Koide, Shibata, 
& Kudoh (1999) and Koide (2003). The detailed presentation of the equation is provided 
in Koide (2003) and will not be repeated here. The numerical scheme is constructed in the 
set of conservation variables D (density), P* (momentum), e (energy), and (for details see 
Koide, Shibata, & Kudoh (1999) and Koide (2003)). The values of the primitive variables v, 
p, Lorentz factor, and p are calculated from P\ e, and B\ Two nonlinear, simultaneous 
algebraic equations with unknown variables, a = 7 — 1 and b = 7(v • B)/c^, are solved using 
a two variable Newton- Raphson iteration method (for details see Koide, Shibata, & Kudoh 

(1999) and Koide (2003)). Here the Lorentz factor is defined 7 = 1/^(1 - Ei=i(Wc)^)- 

In the course of the simulation. Gauss's law for the magnetic field, V • B = 0, and 
the charge conservation law, V^J'* = must be verified. Here, = (cpc, J^, J^, J^) is the 
current four vector. We adapt the Boyer-Lindquist coordinate system (i, r, 9, (p), the metric 
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of Schwarzschild space-time is simplified to 

ho = a, hi = l/a, h2 = r, h3 = rsm9, u!i = (4) 

where the lapse function is a = (1 — rs/r)^/^ and Schwarzschild radius is rs = 2GM/c^ 
(G = 6.67 X 10~^^ Nm^/kg^ is the gravitational constant). Here, r, 6, (j) are radial, latitudinal, 
and azimuthal coordinates, respectively. The line element for the Schwarzschild metric has 
the form 

We use the modified tortoise radial coordinate, r' — ln(r/rs — 1). 

The GRMHD equations described in this section, with the apphed initial and boundary 
conditions (which will be described in the next section), were solved using the simplified 
total variation diminishing (TVD) method. The simplified TVD method was developed by 
Davis (1984) for studying violent phenomena such as shocks (Koide, Nishikawa, & Mutcl 
1996; Nishikawa ct al. 1997, 1998; Koide, Shibata, & Kudoh 1999). This method stems 
from Lax-WendorofF's method, with additional diffusion terms. It was verified in Koide, 
Shibata, & Kudoh (1999) that this method obeys the laws of conservation of energy. The 
conservation laws {dpe/dt + V • a{J + PeCjS) = 0), where a and /3 is lapse function and 
shift vector, respectively, is restraint almost zero within the round-off error with the double 
precision. It was also confrmed that when we take the nonrelativistic hmit {{v/c)^ <^ 1), 
the Newtonian calculation (nonrelativistic MHD simulation) of a magnetized jet resulted 
(Koide, Nishikawa, & Mutel 1996). 

3. 3-D GRMHD Simulations: Simulation Parameters and Numerical 

Techniques 

In this section we present our results showing the fully three dimensional simulations of 
jet formation using general relativistic magnetohydrodynamics. In the assumed initial state 
of our simulations, the region is divided into two parts: a background corona surrounding 
the black hole, and an accretion disk along the equatorial plane of the geometry (see Fig. 
la). The simulation results arc displayed in the Cartesian coordinate transforming from 
the Boyer-Lindquist coordinates (t, r, 9,(j)) coordinate {x = r sin 9 cos (j),y = r sm9 siiKp, z = 
rcos9). The black hole is located at the origin (^x = y — z — 0). The coronal plasma 
is set in a state of transonic free-fall flow (Bondi accretion), with an adiabatic exponent, 
r = 5/3, specific relativistic enthalpy, h — 1.3, and sonic point located at r = 1.6rs. The 
Keplerian disk region is located at r > td = 3rs, |cos^| < S, where S — 1/8, td is the disk 
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radius and rs the Schwarzschild radius. In this region the density is 100 times that of the 
background corona, while the orbital velocity, v^, is relativistic and purely azimuthal: = 
vk = c/[2(r/rs — 1)]^^^, where vk is the Keplerian velocity. (Note that this equation reduces 
to the Newtonian value, = ^jGMjr^ in the non-rclativistic limit rs/r <^ 1). The pressure 
of both the corona and the disk are assumed to be equal to that of the transonic solution. 
The initial conditions for the plasma surrounding the black hole in terms of the transonic 
flow velocity, density and pressure (t'sc, Pffc^Pffc) and the disk properties (fdis, Pdis^Pdis) are: 

P = Pffc + Pdis (6) 

{lOOpffc (r > td and |cot^^| < 5) , . 

(r < re or |cot^^| > S) 

, , r (0,0, Wk) (r > ro and |cot6'| < 5) 

t (-t;ffc,0,0) (r < ro or \coW\>b) 

where we set b = 0.125 and the smoothing length is 0.3rs. The initial conditions for the 
free-falling corona are based on theories by Bondi & Hoyle (1944), Bondi (1952), and Michel 
(1972). 

An initially uniform magnetic field is applied perpendicular to the accretion disk. It is 
set to the Wald solution (Wald 1974), which represents the uniform magnetic field around 
a Schwarzschild black hole; Bj. = BqCosO and Bg = —aBosinO (where a is the lapse function 
defined in section 2). At the inner edge of the accretion disk, the proper Alfven velocity is 
va — 0.015c when Bq — 0.3^/p^. The Alfven velocity in the fiducial observer frame is 

VA = ^ (9) 

^p+[Tp/{T-l) + B^]/c^ 

The plasma beta of the corona, which is defined as the ratio of gas pressure to magnetic 
pressure, is: f3 = p/{B'^/2) = 1.40 at r = 3rs. The sound speed is calculated from vg = 
ciTp/hy/"^. Our simulation domain is, l.lrg < r < 20rg, 0<^^<7r, O<0< 27r, with a 
mesh size of 100 x 60 x 120. The effective linear mesh spacing at r = l.lrg and at r = 20rs 
is 5.38 X 10^'^rs and 0.97rs, respectively. The angular spacing along the polar and azimuthal 
directions is 5.2 x 10~^ rad. At the z-axis {9 — 0, tt) the fixed boundary is used, therefore the 
kink instability is suppressed near the axis. (However, a warping instability near the equator 
is not suppressed.) The following radiative boundary condition is imposed at r = l.lrs and 
at r = 20rs: 

u^+^ = + <+i - <, (10) 

where the superscript n denotes the time step, the subscript refers to the boundary node, 
and the subscript 1 is the node neighboring the boundary. The computations were done on 
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an SGI ORIGIN 2000 system with 0.898 GB of internal memory. 47 hours of CPU time were 
needed for 10,000 time steps. 

Contributions to jet formation from the gas pressure and electromagnetic field compo- 
nents are found using the following equations derived in Koide, Shibata, & Kudoh (1999): 

Wgj> = -V ■ Vp (11) 

Wem = V • (E + J X B) (12) 

These equations will be used for showing the energy sources for generating the jet in Figs. 
6b and 7b. 



4. Simulation Results: Jet Formation 

The evolution of jet formation for the simulation performed in the region previously 
described (l.lrs < r < 20rs, < ^ < tt, and < < 27r with a grid resolution of 
100 X 60 X 120) is shown in Figure 1. The physical values calculated in spherical coordinates 
are transformed to Cartesian coordinates to display their values (for only limited region near 
the black hole (|a;|, \z\ < Srs). The black hole is located at the origin {x = y = z = Q). 
The parameters used in this simulation are the same as those of the axisymmetric simulations 
shown in Fig. 6 of Koidc, Shibata, & Kudoh (1999). Figure 1 shows the evolution of the 
proper density on a logarithmic scale (colored shading) with the flow velocity (shown by the 
directional arrows) in the x — z plane (|/ = 0). The maximum values are normalized for all 
panels. The black circle represents the black hole's Schwarzschild radius. Figure la presents 
the initial conditions described in Section 3. Figure lb {t — 39.2ts) shows that a shock forms 
in the region of a; = 3.5rs. This shock is due to the rapidly infalling material reaching the 
centrifugal barrier located at x = 2rs. At this time, the initial signs of jet formation are 
seen. When the simulation reaches t = 60.0rs (Fig. Ic), the jet structure is clearly seen. The 
formed jet is identified by higher density (Fig. Ic) and pressure (Fig. 2c) along the vertical 
line at a; = ±4rs. The arrows show the velocity of the jet. This phenomenon was also seen 
in previous three-dimensional simulations (Nishikawa at al. 2002, 2003). In the last frame 
(Fig. Id) at t — 128. 9ts, the accretion disk has thickened and the jet appears to fade into a 
weak wind that flows near the accretion disk. Symmetry at the equator is seen as expected 
even with the freedom of the three-dimensional simulation as no asymmetric perturbation 
has been applied. 

The proper pressure, over the same time scale, is shown in Figure 2. The jet forms in 
the shape of a hollow cylinder which is shown clearly in this figure. The pressure in this 
region is high due to the shock and adiabatic heating at the boundary of the decelerated 
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flow. This high pressure contributes greatly to the jet formation (see Equation 11), which 
is shown by the green curves in Figs. 6b and 7b. A bi-conical high pressure region with 
velocity arrows directed away from the disk in Figure 2c represents a jet. However, at the 
later time the jet is quenched and the pressure decreases, as the shock disappears at r = 3rs- 
Figure 2d shows the faded jet and the wind generated along the disk. 

The evolution of the ratio of gas to magnetic pressure (3 — p/{0.5B^) in the x — z 
plane {y = 0) is shown in Fig. 3. In this figure the maximum and minimum values are not 
normalized in order to see the relative values at each time. The initial cold accretion disk 
is heated and the shock is created around x/r^ = ±3 as shown in Fig. 3b. The generated 
conical jet is seen in Fig. 3c. This structure is also found in tori simulations (Fig. 2 in Hirose 
et al. (2004)). The jet is ejected along the twisted magnetic field (shown in Fig. 4c). Two 
layers of high /3 regions are found near the black hole. After the jet is faded, the accretion 
disk is broadened as seen in Fig. 3d. 

The evolution of jet formation is closely related to the deformed magnetic field as shown 
in Fig. 4. The magnetic field lines, which are initially perpendicular to the accretion disk 
(Fig 4a), are twisted by the accretion disk as the simulation progresses and are pinched by 
the falling corona and accretion material (Fig. 4b). Later the disk drags the magnetic field 
further in the azimuthal direction, transferring angular momentum outward even as matter 
falls toward the black hole. The pinching twisting field lines increase the magnetic field 
pressure and generate a jet near the black hole (Fig. 4c). Equation 12 presents the power of 
the electromagnetic forces which drives jet formation as is also shown by the blue curves in 
Figs. 6b and 7b. At t = 128. 9rs the magnetic field has relaxed and is consistent with what 
is seen in Figures Id and 2d where the jet fades. 

Recently, Hirose et al. (2004) have analyzed the magnetic field structures found in 3- 
D GRMHD simulations of tori in the Kerr metric (De Villiers, Hawley, & Krolik 2003). It 
should be noted that their different initial magnetic field conditions lead to different magnetic 
field dynamics. For example, the sample field lines in Fig. 6 in Hirose et al. (2004) for the 
spin parameter a/M — 0.9 (even a/M — 0) are substantially different from the magnetic 
field fines in Fig. 4. Since the magnetized disk in our model is stable to MRI, the magnetic 
fields are never tangled. De Villiers et al. (2004) have found the unbound outfiows in the 
axial funnel region, which consist of two components: a hot, fast, tenuous outfiow in the 
axial funnel proper, and a colder, slower, denser jet along the funnel wall. 

The time required to form the jet in this three-dimensional simulation is similar to the 
two-dimensional axisymmetric simulations. Ultimately, we expect that the dynamics of jet 
formation are modified by the additional coordinate freedom in the azimuthal dimension 
(when axisymmetry with respect to the 2;-axis is removed), such as the difference between 
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axisymmetric 2-D and full 3-D MHD simulations found by Matsumoto (1999). Due to the 
growth of non-axisymmetric instabilities, the magnetic field lines were compressed into helical 
bundles in the jet. In this simulation with a thin disk a global and visible instability was not 
found in the x — y plane. One of the possible reasons is that the strong magnetic field along 
the 2;-axis (low j3 value) prevents growth of MRI. The magnetic fields found with MRl are 
usually along the azimuthal direction (the initial 13 value is larger) (e.g. De Villiers, Hawley, 
& Krolik 2003). 

In order to analyze the properties of jet formation in this simulation, we present the 
physical variables using several onc-dimcnsional slices along the radial axis and the jet for- 
mation axis. Figures Ic and 2c show the density and pressure at time 60. Org corresponding 
to the time of our one-dimensional shoes shown in Figs. 5-7. At this time the generated 
jet has maximum velocity. The parameters for these slices where chosen in order to perform 
comparisons with Koide, Shibata, & Kudoh (1999), and are plotted showing one dimensional 
radial (x-axis) slices along z — z — 5.6rs and x — 4.48rs in the Cartesian coordinate (for 
all figures at y = 0) . 

The parameters along the x-axis {aX, y — z — Q) with the event horizon at x/r^ — ±1.0 
are shown in Figure 5. The central gaps in these panels correspond to the black hole region. 
Figure 5a which shows the density (blue), pressure (green) and ^-component magnetic field 
energy Bl/2 (red) reveals a shock around x = ±3.0rs outside the centrifugal barrier at 
X = ±2.0rs. Shock and adiabatic heating due to strong deceleration of the fiow yields a 
rclativistic region within x = ±2.0rs where p ~ pc^. Figure 5b shows the power contributions 
from the gas pressure (green) and electromagnetic field (blue) components independently as 
depicted in Equations 11 & 12. These energies contribute to decelerating the accretion disk 
plasma (attention must be given to the sign {—W)). The electromagnetic force is produced 
by the magnetic tension caused by the strongly pinched magnetic field lines as shown in Fig. 
4c. The individual magnetic field components are shown in Fig. 5c. The z-component 
(red) is pinched by the falling corona and accretion disk. However, other components {B^ 
(blue) and By (green)) are close to zero. Slight humps in B^ (blue) at a; = ±2rs are found 
near the black hole but further investigation is necessary to ascertain the cause. The three- 
velocity components Vy^ (blue), Vy (green), and (red) are plotted in Fig. 5d. The value of 
Vy reverses at the location of the centrifugal barrier ai x — ±2rs. The inner edge of the disk 
enters into the unstable orbit region, |a;| < 3rs, as shown in Fig. 5d. The material accretes 
to the points around x — ±2rs but is slowed down by the compressed magnetic field {B^ 
(red)) as seen in Fig. 5c. The advection flow in the unstable orbit region stops a,t x ±2rs, 
and the plasma is ejected in the ^-direction as shown in Figs. Ic and 2c. Inside the barrier, 
the plasma is accelerated rapidly to a relativistic velocity near the horizon as seen in Fig. 5d 
(vx (blue)). This highly heated, high density region should contributes to radiation which 
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can be observed by present and planned space based observatories (see Fabian et al. 2000; 
Reynolds and Nowak 2003). 

Figure 6 shows the same parameters as in Figure 5, but along the z — 5.6rs slice 
(y = 0) . Since this shce is well above the Schwarzschild radius of the black hole, there is not 
a central gap in the figures. Figure 6b clearly shows the separate layers of the jet resulting 
from the gas pressure (green around x = it4rs) and the electromagnetic forces (blue around 
X = ±5rs). A double-layered jet is found where the inner part of the jet is a gas-pressure- 
driven jet and the outer one is a magnetically driven jet (Koide, Shibata, & Kudoh 1999). 
From the velocity profiles in Figure 6d, the jet appears to have a width from 3rs < |.x| < 5rs. 
The electromaguctically driven component has the larger magnetic field By, and a large 
azimuthal velocity Vy. It should be noted that at this time (ts = 60.0) a torsional Alfven 
wave is excited (Koide, Shibata, & Kudoh 1999). Later the outwards propagation of the 
TAW contributes to formation of a wind, and as will be discussed later, the peak of will 
move outward. 

Figure 7 presents the physical variables along the x — 4.48rs shce (|/ = 0). This slice 
is along the vertical axis of the jet rather than the horizontal plane of the disk as in the 
previous figures (Figs. 5 and 6). The density (blue) and pressure (green) increases are 
expected around the accretion disk in Fig. 7a. The equatorial plane is located at 2; = and 
the density in the disk is much higher than other values as shown in Fig. 7a. Figure 7b shows 
the contribution to jet generation due to the electromagnetic field (2 < |-2/rs| < 4) (outer 
part) and gas pressure (5 < \z/rQ,\ < 8) (inner part), and is consistent with the increase 
of the jet Vz (red) in Fig. 7c. However, the region of acceleration due to gas pressure is 
smaller than in the axisymmetric case. This may come from the fact that the time here is 
somewhat later than the t — 52ts in the axisymmetric case and/or due to 3-D effects. The 
three- velocity components are plotted in Fig. 7c. It should be noted that the jet has a larger 
azimuthal component (vy) near the disk than the z component (f^) which becomes larger 
far from the disk. A slightly positive radial velocity fx (blue) is found far from the disk in 
Fig. 7c, and this indicates outward angular momentum transfer. The Alfven velocity va and 
the sound velocity fg are shown along with the absolute value of in Fig. 7d. The velocity 
Vz gradually increases in the magnetically driven region 5rs < z < 7rs but does not exceed 
either the sound or Alfven velocities. 

Figure 8 shows one dimensional slices of the velocities (blue), Vy (green), and Vz 
(red), at ^ = 5.6rs (y = 0) for times (a) t = COrg, (b) t = 39.2rs, (c) t = 60.0rg, and (d) 
t = 128. 9rs. The jet is seen in Vz (red curve) and its maximum velocity can be found in Fig. 
8c. Comparison between the velocity structure at t = 60ts (Fig. 8c) and 80rs (not shown) 
reveals that the jet velocity {vz) and the azimuthal velocity {vy) are smaller at t — 80ts. 
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Additionally, the location of peak velocities (v^) are shifted outward at the later time. This 
may come from the fact that the wind is caused by the unwinding magnetic fields with the 
TAWs which are propagating outward in the disk. These results are consistent with the fact 
that the jet has faded and a weak wind is generated around x = 7rs as seen in Fig. 8d. 
Furthermore, the rotation of the accretion disk (vy) becomes very small and rotates slightly 
in the reverse direction around x = ±3.5rs, and ±7rs. We plan to investigate how these 
changes are affected by inflow of matter to the system in the future. 

A comparison of our present results with previous two-dimensional results shown in 
Koidc, Shibata, & Kudoh (1999) reveals some minor differences that may be due to the 
three-dimensional modeling. The most prominent difference between this 3-D simulation 
and previous 2-D simulations lies in the longer run time and hence shows a significant 
change in the disk-jet system to a disk- wind system (Fender 2002, 2003). This transition 
may be caused by outward TAW propagation in the disk. The jet is mainly produced by 
the gradient of gas pressure enhanced by the shock around r = 3rs, as shown in Figs. Ic, 
3c and 6b. After t = 60rs, the TAW moves outward in the disk and is weakened. The wind 
appears while the TAW is moving outward. After outward propagation of the TAW, the 
pressure of the outer disk becomes larger than the initial condition (Fig. 2a). Therefore, 
the wind may be caused by unwinding of the twisted and pinched magnetic fields due to the 
decreased accreting matter, which creates the TAW propagating outward in the disk. This 
transient change may be controlled by the mass accretion rate onto the black hole. It should 
be noted that in the simulations with initial weak azimuthal magnetic fields the magnetic 
fields are tangled as shown in Fig. 6 in Hirose et al. (2004) as a result of MRI where TAWs 
are not excited. 



5. Summary and Discussion 

We have presented a fully 3-D GRMHD simulation of a Schwarzschild black hole accre- 
tion disk system with moderate resolution. Because no deliberate asymmetric perturbation 
was introduced, the present 3D simulation remains largely axisymmetric. While we found 
some numerical differences between this simulation and previous two-dimensional simula- 
tions, our results reaffirm the previous result of a two-layered jet. Our results suggest that 
the jet develops as in the 2-D case. Our longer simulation has allowed us to study jet evolu- 
tion. At the end of this simulation the jet fades and a weak wind is generated by a thickened 
accretion disk. This phenomena was not observed in the two-dimensional simulations be- 
cause they did not run for this duration. 

In this simulation at the outer boundary of the accretion disk no matter is injected. 
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therefore after accreted matter is ejected from the accretion as a jet, the power to generate a 
jet is dissipated and the jet is switched to a wind. However, if streaming matter is injected at 
the outer boundary of accretion disk, transient changes may be controlled by accretion rates 
with instabilities and be related to the state transitions. Such disk-jet coupling in black hole 
binaries is reviewed by Fender (2002, 2003) and Fender, Belloni, & Gallo (2004). Black hole 
binaries exhibit several different kinds of X-ray 'state'. The two most diametrically opposed, 
which illustrate the relation of jet formation to accretion, can be classified as low/hard/off 
and high/soft states (Fender 2003). These two states provide different luminosity. Pellegrini 
et al. (2003) have discussed a nuclear bolometric luminosity and an accretion luminosity Lace 
in terms of the accretion rate M and jet power. Clearly these issues can be investigated by 
further 3-D GRMHD simulations, and future simulations will investigate jet formation with 
different states of the black hole including streaming material from an accompanying star 
(e.g.. Miller et al. 2003). 

Additionally, we plan to study the development of instabilities along the azimuthal 
direction such as the proposed accretion-ejection instability (AEI) (e.g. Tagger & Pellat 
1999; Gaunt & Tagger 2001; Rodriguez et al. 2002; Vaniere, Rodriguez & Tagger 2002), the 
magnet orot at ional instability (MRl) (e.g. Balbus & Hawley 1991, 1998), and the "screw" 
(current- driven helical kink) instability (e.g. Li 2000). Such simulation studies may indicate 
how the observed variabihties in AGN jets (Blandford 2003; Mirabel 2003) and microquasars 
(e.g. Eikenberry et al. 1998; Ueda et al. 2002; Nobil 2003; Rau, Greiner, & McGoUough 
2003) can be produced. Initial magnetic field configurations affect the dynamics of accretion 
disks and the associated phenomena such as jet formation and its variabilities, and thus we 
plan to study different initial field configurations. 

The Fe Ka fluorescent emission hne in active galactic nuclei (AGN) is beheved to origi- 
nate in a relativistic accretion disk around a black hole (e.g. see reviews Fabian et al. 2000; 
Reynolds & Nowak 2003). The X-ray continuum emission of MGG-6-30-15 is highly variable 
(e.g. Vaughan & Fabian 2003). The source spectrum and variablity can be explained by 
two- components consisting of a variable power law component (PLC) and an almost con- 
stant reflection-dominated component (RDC) containing the iron line (Fabian & Vaughan 
2003). The region where the flare (or multiple flares) which created the PLC must be located 
close to the rotation axis and within 3-4 gravitational radii to reproduce the highly variable 
behavior, (Miniutti et al. 2003). Since this behavior is most likely powered by magnetic 
field from the disk, possibly linking to the hole (e.g. Blandford & Znajek 1977; Wilms et al. 
2001), and perhaps forming the base of jet, it is essential to perform GRMHD simulations. 
These issues are being investigated with 3-D GRMHD simulations with a Kerr black hole 
with streaming matter from an accompanying star. 
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Fig. 1. — The panels present the time evolution of the proper mass density (color) with 
velocity arrows {vx,Vz) in a transonic free-fall (steady-state falling) corona with an initially 
uniform magnetic field, at (a) t = O.Org, (b) t = 39.2rs, (c) t = 60.0rs, and (d) t = 128. Org. 
The black hole's Schwarzschild radius is presented by the central black circle in each panel. 
The jet is fully formed at t = 60. Org (c). At the later time the wind is formed with a wider 
angle. 
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Fig. 2. — These panels present the time evolution of the proper pressure (color) with velocity 
arrows {v^, Vz) in a transonic free-fall (steady-state falling) corona with an initially uniform 
magnetic field, at (a) t = O.Ors, (b) t = 39.2rs, (c) t = 60.0rs, and (d) t = 128.9rs. 




Fig. 3. — The time evolution of gas to magnetic pressure (3 = p/ (0.5-B^) (color) with velocity 
arrows at (a) t = O.Ors, (b) t = 39.2rs, (c) t = 60.0rg, and (d) t = 128.9rs. The profiles of 
generated jets are seen clearly in Fig. 3c. 



Fig. 4. — These panels show the time evolution of the magnetic field lines associated with a 
Schwarzschild black hole and an accretion disk. The field lines are traced from the bottom 
{z = — 8rs) with the different radius (red: r = 2.67rs; blue: r = 5.33rs) at (a) t = O.Ors, (b) 
t = 39.2rs, (c) t = eO.Ors, and (d) t = 128.9rs. 




Fig. 5. — This plot represents the one dimensional slice through the simulation domain at 
z = a.t t = 60.0rs. The location of the black hole is excluded in the plots. Frame (a) shows 
the proper mass density p (blue), proper pressure p (green), and magnetic field energy 
(red). Frame (b) shows the power —Wgp (green) and —Wem (blue) associated with the gas 
pressure force and electromagnetic force which contribute to decelerating the accretion disk 
plasma. Frame (c) shows the individual magnetic field components (red) and Bj (green), 
and Sx (blue). Frame (d) shows the velocity components fx (blue), Vz (red), and Vy (green). 




Fig. 6. — This plot represents various physical quantities on the z = 5.6rs surface at t = 
60.0rs. Frame (a) shows the values for proper mass density p (blue), proper pressure p (green) 
and z component magnetic field energy (red). Frame (b) is the power contribution from 

the pressure, Wgp (green) and the electromagnetic force, Wem (blue). We can find the two- 
layer acceleration region in the jet. Frame (c) shows the individual magnetic field components 
i?x (blue). By (green), and (red). Frame (d) has the components of velocity (blue), Vy 
(green), and Vz (red). 




Fig. 7. — This plot represents various physical quantities on the x = 4.48rg surface at 
t = 60.0rs. Frame (a) shows the values for proper mass density p (blue), proper pressure 
p (green), z component magnetic field energy (red), and perpendicular magnetic field 

energy + (yellow). Frame (b) is the power contribution of the pressure, Wgp (green) 
and the electromagnetic force, Wem (blue). Frame (c) shows the components of the velocity 
Vy (green), (red), and (blue). Frame (d) shows the Alfven velocity va (blue), the sound 
velocity fs (green), and the velocity \vz\ (red). 




Fig. 8. — The time evolution of one dimensional slice of velocities through the simulation 
domain at z = 5.6rs is plotted at (a) t = O.Ors, (b) t = 39.2rs, (c) t = 60.0rs, and (d) 
t = 128. 9rs. The blue, green, and red curves present Vx,Vy, and v^, respectively. 



